{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import numpy as np\r\n",
    "import pandas as pd\r\n",
    "import statsmodels.api as sm\r\n",
    "import seaborn as sns\r\n",
    "from pymongo import MongoClient"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "# 输入是一DataFrame，每一列是一支股票在每一日的价格\r\n",
    "def find_cointegrated_pairs(dataframe):\r\n",
    "    # 得到DataFrame长度\r\n",
    "    n = dataframe.shape[1]\r\n",
    "    # 初始化p值矩阵\r\n",
    "    pvalue_matrix = np.ones((n, n))\r\n",
    "    # 抽取列的名称\r\n",
    "    keys = dataframe.keys()\r\n",
    "    # 初始化强协整组\r\n",
    "    pairs = []\r\n",
    "    # 对于每一个i\r\n",
    "    for i in range(n):\r\n",
    "        # 对于大于i的j\r\n",
    "        for j in range(i+1, n):\r\n",
    "            # 获取相应的两只股票的价格Series\r\n",
    "            stock1 = dataframe[keys[i]]\r\n",
    "            stock2 = dataframe[keys[j]]\r\n",
    "            # 分析它们的协整关系\r\n",
    "            result = sm.tsa.stattools.coint(stock1, stock2)\r\n",
    "            # 取出并记录p值\r\n",
    "            pvalue = result[1]\r\n",
    "            pvalue_matrix[i, j] = pvalue\r\n",
    "            # 如果p值小于0.05\r\n",
    "            if pvalue < 0.05:\r\n",
    "                # 记录股票对和相应的p值\r\n",
    "                pairs.append((keys[i], keys[j], pvalue))\r\n",
    "    # 返回结果\r\n",
    "    return pvalue_matrix, pairs"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "from datetime import datetime\r\n",
    "# 时间处理\r\n",
    "def loadDate(str):\r\n",
    "    date = datetime.strptime(str, '%Y%m%d').strftime('%Y-%m-%d')\r\n",
    "    return (date)\r\n",
    "# 数据查询\r\n",
    "def selectSotck(adjustflag,codeList, startDate, endDate):\r\n",
    "    client = MongoClient('47.100.19.231', 27018)\r\n",
    "    db = client.admin\r\n",
    "    db.authenticate(\"quant\", \"Qweasd123\")\r\n",
    "    db = client.quant\r\n",
    "    mycol = db.testKline\r\n",
    "    \r\n",
    "    data_list = []\r\n",
    "    for item in codeList:\r\n",
    "        myquery = {\r\n",
    "            \"adjustflag\": adjustflag,  # 复权类型\r\n",
    "            \"code\": item,  # 股票代码\r\n",
    "            \"date\": {'$gt': startDate, '$lt': endDate}  # 时间区间\r\n",
    "        }\r\n",
    "        mydoc = mycol.find(myquery)\r\n",
    "        for x in mydoc:\r\n",
    "            data_list.append(x)\r\n",
    "    fields = ['date','code','close']\r\n",
    "    df = pd.DataFrame(data_list, columns=fields)\r\n",
    "    df.columns = ['date','instrument','close']\r\n",
    "    df['date'] = df['date'].map(loadDate)\r\n",
    "    df['close'] = df['close'].astype('float64')\r\n",
    "    return df"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "instruments = [\"sz.002142\", \"sh.600000\", \"sh.600015\", \"sh.600016\", \"sh.600036\", \"sh.601009\",\r\n",
    "              \"sh.601166\", \"sh.601169\", \"sh.601328\", \"sh.601398\", \"sh.601988\", \"sh.601998\"]\r\n",
    "\r\n",
    "# 确定起始时间\r\n",
    "start_date = '20150101' \r\n",
    "# 确定结束时间\r\n",
    "end_date = '20170719' \r\n",
    "# 获取股票总市值数据，返回DataFrame数据格式\r\n",
    "prices_temp = selectSotck('1',instruments, start_date,end_date)\r\n",
    "\r\n",
    "prices_df=pd.pivot_table(prices_temp, values='close', index=['date'], columns=['instrument'])\r\n",
    "pvalues, pairs = find_cointegrated_pairs(prices_df)\r\n",
    "# #画协整检验热度图，输出pvalue < 0.05的股票对\r\n",
    "sns.heatmap(1-pvalues, xticklabels=instruments, yticklabels=instruments, cmap='RdYlGn_r', mask = (pvalues == 1))\r\n",
    "print(pairs)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[('sh.601009', 'sh.601169', 0.02365346158056714), ('sh.601328', 'sh.601988', 0.015970238759557664), ('sh.601328', 'sh.601998', 0.021156769414874762), ('sh.601988', 'sh.601998', 0.0474966733666659)]\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEjCAYAAAAomJYLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debwcVZ338c83CTsii6CoYEDZVyEw6IOCogziDIsOIjgjLpjBdYRxwUEdHXDYdGZ0RAIiDvOoIMMue8wg4KOBsISEEAibQtzYQUSW5P6eP865pHPpe2+lu0+6qvm+edUrfauqv326E/rcqlN1fooIzMzMRprQ7waYmVk9uYMwM7O23EGYmVlb7iDMzKwtdxBmZtbWpH43oGaKXNL1s98cVSKWbdbZukguwNqPPlokd/o23yySu+fvv1UkF+DO/3NEkdw1rvpskdy1V16/SC7A0L+fUiR3pcM/XiQ3Ft5SJBdAGx6hbjN+pM0qf+ccHHd0/XrLykcQZmbWlo8gzMz6ZFLNv4Fr3jwzs8E1oebncNxBmJn1iTsIMzNra8JyH3ZeNsX6LyXfknSXpDmSdmjZtpekO/K2I1vWnyjp9rz/+ZLWzOvXkXSVpCclfXuU17tI0q0tPx8h6bacNUPSa0q9VzOzTkyYUH3pS/sKZr8D2CQvU4GTASRNBE7K27cEDpK0ZX7OdGDriNgWWAB8Ia9/GvgS8Jl2LyTpXcCTI1bfDEzJWecAJ/TmbZmZ9cZAdBCSVpN0iaRbJN0q6fOSZudlrqR21/LuC/x3JDOBNSWtD+wM3BUR90TEs8BZeV8i4sqIWJSfPxN4dV7/p4j4OamjGNm21YEjgGNa10fEVRHx1MgsM7O6mDSp+tIPVfulvYDfRsR2EbE1MC0ito+I7YHLga+3ec6rgPtbfl6Y1422fqQPAZdVaNvRwDeAp8bY58OjZUmaKukGSTeceuqpFV7OzKw36n4EUbVfmgt8XdLxwMURcS2ApPcAOwB7tnlOu+GXGGP9kidKRwGLgB+O1ShJ2wOvi4jDJU0eZZ+/BaYAu7XbHhGnAsM9g4tjmNlyMxBXMUXEAkk7AnsDx0q6EjgX+Crw5ohY3OZpC4ENWn5+NfBbYMVR1gMg6RDgr4A9YvxqRm8AdpT0q/xe1pP0s4jYPWe9DTgK2C0inqnyXs3Mlpe6dxBVxyBeCTwVET8gnU56C2ns4P0R8eAoT7sIeH++mmkX4PGI+B0wC9hE0kaSVgTem/dF0l7A54F9WsYPRhURJ0fEKyNiMrArsKClc3g9cErOeqDK+zQzW54kVV76oeoppm2AEyUNAc8BFwP/CHx3uOERsb2kw/LjacClpCOOu0jjAx/M2xZJ+gRwBTAROD0i5uXX+TawEjA9586MiMMA8lHCGsCKkvYD9oyI28Zo84nA6sD/5Kz7ImKfiu/XzKy4gZhqIyKuIH2ht/pqm/2mtTwOoO0UjRFxKakDGbn+dWO0YfI4bfwVsHXLz28ba38zs36r+ymmmvdfZmaDyx2EmZm15Q7CzMzacgdhxXz9pquKZf/rWhsWyd3zd/9RJPfy9T5VJBdgl98dM/5OHZi/+ZeK5E584zpFcgFWfdVLiuTeu0eZqn2Tv/3uIrm94g7C2P1VXyuSe8WvpxbJNbPlY9LEfrdgbO4gzMz6xEcQZmbWljsIMzNryx2EmZm1NaHmJeU66r8k/UrSyyrst6GkKyXNz9XdJuf1G0m6TtKdkn6c52TqtArd2pKm56zpktZq2faFvP8dkv6yk/dqZlaKJqry0g+lD3D+GzgxIrYgFQoanjTveODfI2IT4FFSvQborArdkcCMnDUj/0ze/l5gK1I9i+/kHDOzWpgwaULlpS/tG2+HNtXkDsybPinpplxRbvM2z9sSmBQR0wEi4smIeEpp5ry3ksqAApwB7JcfL3MVuvznGaNknRURz0TEvaRJA3eu+LmYmRU3CEcQI6vJXZ7XPxQRO5B+y29XK3pT4DFJ50m6WdKJ+Tf4dYDHWkqLtlaU66QK3cvzNOLkP9cbJ2sprihnZv2iCaq89EOVQeoXVJPL02efl7ffCLxrlOw3Aa8H7gN+DHyAXPthhOHCQB1XoWuj0nNcUc7M+qVfRwZVjXsEERELgB1JHcWxkr6cNw1XaFtM+45mIXBzPi20CLiAVJ70IdKpo+HntFaUG60K3WjrAf6QT0OR/xwe5xjrOWZmfTdhgiovfWnfeDu0qSa3wzhPGTYLWEvSuvnntwK35ToRVwF/k9cfAlyYHy9zFbr85yGjZL1X0kqSNiINfF9fse1mZsVNWGFC5aUfqpxiGllN7qMsGWBeiqQpwGERcWhELJb0GWBGHpi+Efhu3vXzwFmSjgFuBr6X13dShe444GxJHyadyjogP2eepLOB24BFwMdHqZ1tZtYXqvmdcuN2EKNUk5vcsv0GYPeWx4e2bJsObNsm8x7aXFHUYRW6h4E9RnnO14AyM+WZmXWp7mMQvpPazKxP+nV1UlXuIMzM+sQdhJmZteVTTNZIz5w3q0juilvdWyR38t2Hjr9Th9a47PLxd+rALlceMv5OnVhz3fH36VC8bKMiuRs9+UiR3KHrrimSC6Ctus+Y2KcpNKpyB9Fgx76x3J3fT1+17/g7mVlX6n4EUe/uy8xsgPVyqo3RZrxu2f5SST/J8+rNk/TB8TJ9BGFm1ie9OoJomfH67aRZJGZJuigibmvZ7eOkm5X/Ot/AfIekH+YJUNtyB2Fm1ic9vFHu+RmvASQNz3jd2kEE8JJ84/LqwCOkm4hH5VNMZmZ9siynmFpnns7L1JaoKrNXfxvYgjQn3VzgHyJiaKz2DUJFuQPy+bShPNXH8PrJkv4saXZepnXyXs3MSpm4woTKS0ScGhFTWpbWq1SqzF79l8Bs4JXA9sC3Ja0xVvsGoaLcraTpxttdz3Z3RGyfl8N6/ebMzLrRw4JBVWav/iBwXi7IdhdwL/CCYm+tGl9RLiLmR8Qd470PM7O66eFVTGPNeD3sPvK8dZJeDmwG3DNW6CBUlBvLRvm1r5b0pgr7m5ktPxNVfRlD/j4dnvF6PnB2ntH6MEnDZ0+OBt4oaS4wA/h8RDw0Vu4gV5T7HbBhRDwsaUfgAklbRcQTrTvlgZ6pAKeccgpTp05tE2Vm1nu9nIup3YzXETGt5fFvgT2XJbPKdN8L8hfs3qSKclfmTZUrygFIugDYBTidXFEu93pVKsqtOMr6sdr9zHAbI+JGSXeTjmpuGLGfS46aWX+sMLHfLRjTIFSUG63d6+ZTWkjamDTwPeb5NjOz5amHg9RFVBmD2Aa4XtJs4CjgmNF2lDRF0mkAuXrbcEW5uaTTRK0V5Y6QdBdpTKK1otw9pIpy3wU+lrPanl/Lr7m/pIXAG4BLJA0XN3ozMEfSLaQB8cMiosyMYGZmnZig6ksfDEJFufOB89usPxc4t12WmVkt1HyyPk+1YWbWJy4YZGZm7U2s92xH7iDMzPpEK7iDMHvepOm3jb9TB14y9/7xd+rQE799YvydOnHFvCKx53/2DUVyAd71le+Nv1ON7PeOTYtlnzf+LuPzGIQ10cpHXTj+Tp34eLkvL7Om8RiEmZm15yMIMzNrq3cFg4pwB2Fm1icepDYzs/Z8isnMzNqp+yD1wJYczdu2lfTLvH2upJU7eb9mZkVMnFB96YOBLTkqaRLwA9IkfVuR5ot6rvdv0cysQz0qGFTKIJcc3ROYExG35P0ezjPMmpnVQg9LjhYxyCVHNwVC0hW5I/tcu50kTZV0g6QbTj311Ha7mJmVscKE6ksfDHLJ0UnArsBOwFOkuhQ3RsSMpUJcUc7M+qRfhYCqGrdbiogFwI6kjuJYSV/OmyqXHM1HCxeQqtE9RC45mverUnJ0tPVjWQhcHREPRcRTpFoSVavhmZmVN2FC9aUfzRtvh6aWHCUVOdpW0qq5M9oNKDNTnJlZJ2peUW5gS45GxKPAv5E6l9nATRFxSYX3a2a2fNT8CGJgS47mbT8gXepqZlY/novJzMzamjSx3y0YkzsIM7N+8RGEWXmP7ffGYtn/s8OZRXJXXvDOIrnvvP7hIrkAd/7rgePv1IEH//xYkdyHnv5jkdyeqflcTO4gbLkaOumXRXJj7leK5JoV5SMIMzNryx2EmZm15UFqMzNry0cQZmbWTp7XrrbcQZiZ9UvNjyAGoaLc0Xnf2fm1XtmyzRXlzKy+aj7VxiBUlDsxIraNiO2Bi4Ev5+e4opyZ1VvTJ+trQEW5J1pedjWW1HRwRTkzq7dJE6svfTAQFeUkfU3S/cD7yEcQuKKcmdXdAJximgu8TdLxkt4UEY/n9a0V5Sa3ed5wRbnPkKq6bUyqKDdWdbiOKspFxFERsQHwQ9K04MOvvyup09gV2F/SHi8IiTg1IqZExJSpU6e2eRkzs0J62EGMNk47Yp/d83jtPElXj9u88XZoWEW5HwHvbslyRTkzq68ejUGMM047vM+awHeAffK47AHjNm+8HepeUU7SJi2vuQ9we37sinJmVm+9O4IYdZy2xcHAeRFxH0BEPMA4Gl9RDjguD57PIQ1M/0N+jivKmVm9LUMH0TpempfWc+JjjtNmm5J+af+ZpBslvX+85g1CRbl3t9l9eJsryplZfS3D1UkRcSow2pU0Y47TDr8aabhgD2AV4JeSZuZhhPbNq9w6MzPrLfXs6qQq47QLSVef/gn4k6RrgO2AUTuIet/nbWY2yDSh+jK2UcdpW1wIvEnSJEmrAn9BOmU/Kh9BmI1j9/3XLpK7Lq8oknvPf04vkgvwqjseKZK77m+eLJL7h9njjsN27rpvdJ/Ro/sbImKRpOFx2onA6RExT9Jhefu0iJgv6XJgDjAEnBYRt46V6w7CBoK2+UrB9JFDcGY90rtTTG3HaSNi2oifTwROrJrpDsLMrF8m1PsruN6tMzMbZDWf7tsdhJlZv/TwFFMJ7iDMzPrFHYSZmbVV8w6ibhXlNs8V4J6R9JkRWaNVlFtb0vScNV3SWnn9ipK+n+tV3CJp907eq5lZMb27D6KIulWUewT4FGlSwOeNM1PhkcCMnDUj/wzwEYCI2AZ4O/ANqebdtZm9qGjiCpWXfqhVRbmIeCAiZvHC0qBjzVS4b85YKovUkcwYzgUeA6aM937NzJabATiCWJ4V5UYz1kyFL89TgpP/XC+vvwXYN99WvhFpkqrWuUoAV5Qzsz6qeQdRZZB6LvB1SccDF0fEtekgYKmKcu8aJftNwOuB+4AfkyrKjZwfBF446+BIVWYqHOl0YAvgBuDXwC+ARSN3GjFD4niZZma9U/Oz3lWm+14gaUdgb1JFuSvzpsoV5QAkXQDsQvriXlPSpHwUMVp1uJFZo81U+AdJ60fE7yStTx7nyNmHDz9B0i+AO8d7v2Zmy03Nb5SrW0W5sbJGm6nwopyxVFauJLdafvx2YFFEuKKcmdXHAJxi2gY4UdIQafD4oywZYF6KpCnAYRFxaEQszpeqzsgD0zeydEW5syQdA9xMrign6RWkU0JrAEOSPg1sGRFPtJupMGcdB5wt6cOkU1nDdVbXA67I7f4N8HfVPhIzs+Wk6XMxLeeKcr8nnT5q147RKso9TKqQNHL9r4DN2mWZmdVCzU8x1bv7MjMbZE0fpDYzs0LcQZhZO8c8eHeR3AMvPmT8nTr08z8+XCT37scfK5L7/nNWK5LbM+4gzJptvbN/WSb4+r8vk2uNEcswBtHuZrDS3EGYmfVJxFDlfdWHHsIdhJlZnwwtQwcxwR2EmdmLR1C9g+gHdxBmZn2yLEcQ/eAOwsysT+p+BNGkinKnS3pA0q1tXueTudrcPEkntKzfNufNy3UrVu7k/ZqZlbB4aFHlpR8aUVEu+y9SbYqlSHoLqWjQthGx1fBzJU0CfkCaG2or0nQgIwsRmZn1zdAy/NcPTakoR0RcQ+pARvoocFxEPDOckdfvCcyJiFvy+ocjYvF479fMbHmJGKq89ENTKsqNZVPgTfmU1dWSdmpZH5KuyB3Z59o92RXlzKxfhmKo8tIPTakoN5ZJwFqkYkQ7kab+3jiv3zWve4o07fiNETFjqRd2RTkz65PGD1JHxAJSPee5pIpyX86bKleUy0cLF5CKDT1EriiX96tSUW4sC4HzIrkeGAJeltdfHREPRcRTpKnCqxY7MjMrbnE8V3nph6ZUlBvLBTkbSZsCK5I6oSuAbXNluUnAboAryplZbdT9FFOVMYhtgOslzQaOAo4ZbUdJUySdBpAHhIcrys0lzTXVWlHuCEl3kcYknq8oJ2khcATwRUkLJa2Rt50J/BLYLK8fvvLpdGDjfPnrWcAh+WjiUeDfSB3VbOCmiLik2sdiZlZe3Qepm1RR7qBR1j8L/O0o235AutTVzKx2+nX5alW+k9rMrE/6dWRQlTsIM7M+8VxMZmbW1uLozxQaVbmDMBswm/zhmfF36tANm583/k4deGmhSX9WefKE8Xfqo17eByFpL+CbwETgtIg4bpT9dgJmAgdGxDnt9hnmDsKsT7608ylFcuOOrxXJtd7r1SmmPEvFScDbSfeAzZJ0UUTc1ma/43nhhUdt1btitpnZAOvhZa47A3flG5OfJV3yv2+b/T4JnMuSiVPH5A7CzKxPlmU219Z54/IytSXqVcD9LT+/YI47Sa8C9gemVW2fTzGZmfXJslzmOmLeuJHaVaweObfcfwCfj4jFeT69cbmDMDPrk0VDPatAsBDYoOXndnPcTQHOyp3Dy4C9JS2KiAtGCx3YinKSVpT0/Vyv4hZJu3fyXs3MShmKqLyMYxawSf5uXRF4LyNmzo6IjSJickRMJtXj+dhYnQMMcEU54CMAEbENaWT/G5I85mJmtTFEVF7GkmfM/gTp6qT5wNkRMU/SYZIO67R9455ikrQacDbpkGUicHTe9ElJfw2sABwQEbePeN4LKsrl9cMV5Q7Ou54BfAU4OVeDe0DSO0e2IyKuGT4CGWG0inJbAjOG10l6jHSIdf1479nMbHmocGRQWURcSipr0Lqu7YB0RHygSuYgV5S7BdhX0iRJG5FqWmww8smuKGdm/VL36b4HuaLc6cAWwA3Ar4FfAC+4r90V5cysXxYNNXwupohYIGlHYG9SRbkr86bKFeUAJF1A+hI/nVxRLh9F9KyiHKluxRDwsoh4EDh8eCdJvwDu7OJ1zMx6qu6T9Q1sRblcSW61vP7twKKRt52bmfVTD69iKmJgK8oB6wE3SZqfX+/vqn0kZmbLR6+uYiplYCvKRcSvgM3aPcfMrA7qforJd1KbmfVJv04dVeUOwsysTxp/FZOZmZXhU0xmNjDed9haRXInrVzmq+jAX80tkgugLdqVW1g2PsVkZsuVNjuqYPpPCma/+LiDMDOztoZ6WJO6BHcQZmZ94iMIMzNry1cxmZlZW3W/iqkRFeUkrSzp+lwZbp6kr7ZsO1HS7ZLmSDpf0pp5/QqSzsgV5eZL+kIn79XMrJRBmIupG72qKPcM8NaI2A7YHthL0i5523Rg64jYFlgADHcEBwAr5YpyOwJ/P0rBITOzvmh8ByFpNUmX5N/eb5V0YN70SUk35d/QN2/zvBdUlIuIp1oqyp2Tdz0D2C/v80BEzAKea82K5Mn84wp5ibztypbiQzNZMpdTAKtJmgSsAjwLPDHuJ2Jmtpw0voOgJhXlJE3MM8o+AEyPiOva7PYh4LL8+BzgT8DvSAWLvh4Rj7TJdUU5M+uL54ai8tIPjakol6cP3z6PMZwvaeuIuHV4u6SjSBXjfphX7UwqZvRKUsW5ayX9dLiAUUuuK8qZWV8srvk3zrhHEBGxgHQOfy6potyX86bKFeXy0cIFpGJDD5EryuX9lqmiXEQ8BvyMdGQDgKRDgL8C3pdrQQAcDFweEc9FxAPA/wOmVH0dM7PShqL60g+NqCgnad2Wq5NWAd4G3J5/3otUEGifiHiq5Wn3AW9Vshqp3OntFdtuZlbc4ojKSz80paLc+sBVkuaQOp7pEXFxzvo28BJguqTZkqbl9ScBqwO35ud8PyLmVHi/ZmbLRd2PIJpSUW4OaSyjXfteN8r6J0mXupqZ1VLdxyB8J7WZWZ/06+qkqtxBmJn1Sc2nYnIHYWbWLz7FZGY2jgu//2CR3Nj3oSK5ANqi+4yan2FyB2Fm1cXJM4vkXnrmZkVy665fl69W5Q7CzKxPfARhZmZt+SomMzNry4PUZmbWVs0PIAaiotx2+TlzJf0kT83hinJmVnuLo/rSD4NQUe404MhcOe584LN5vSvKmVmtDQ1F5aUfGl9RDtgMuCY/ng68e/hpuKKcmdXYc0PVl34YhIpytwL75McHABvkx64oZ2a11stTTJL2knSHpLskHdlm+/skzcnLLyRtN15mlQ5iLvA2ScdLelNEPJ7Xt1aUm9zmecMV5T4D7ARsTKoopzb7VqooFxHbk2Z73VnS1nnTh4CPS7qRNO33s3l9a0W5jYB/lLRxm9xTI2JKREyZOnXqeM0wM+uZXtWkzr98nwS8A9gSOCifxWl1L7BbRGwLHM2SSpqjanxFuYi4PSL2jIgdgTOBu/OurihnZrXWwyOInYG78vfts8BZwL6tO0TELyLi0fzjTNqXVljKIFSUWy//OQH4IjBcMMgV5cys1palolzr6fC8tJ7yeBVwf8vP4526/zBw2Xjtq3IfxDbAiZKGSIPHH2XJAPNSJE0BDouIQyNicb5UdUYemL6RpSvKnSXpGOBmWirKATcAawBDkj5NOlxaHzgjH0ZNAM5uqSh3kKSP58fnAd/Pj0/Kj28lndZyRTkzq5XFyzD4HBGnMvppocqn7iW9hdRB7Dreaw5CRblvAt9ss94V5cys1p5dlh5ibAtZcoEOjHLqXtK2pFsD3hERD48X6jupzcz6pIezuc4CNpG0EfAb4L2kcdjnSdqQdJbl7/LY8rjcQZiZ9UmvDiAiYpGkT5DO9kwETo+IeZIOy9unAV8m3WbwnXTWn0URMeaFO+4gzMz6pJf1ICLiUuDSEeumtTw+lJYhgCrcQZjZwHpmxp3FslfZvfuMxTWfrc8dhJn13d6P3VEk989femeR3F55tubzfbuDMDPrE5ccNTOzttxBmJlZWx6DMDOztmo+BNGMinJ52z/kehTz8hQcw+u3lzRT0uw8P8nOeb0ryplZrS0eispLPzSiolye2vsjOWM74K8kbZI3nwB8NU8F/uX8M7iinJnV3LOLhyov/dCIinLAFsDMiHgqTx1+NbB/3hakyf0AXsqS+UdcUc7Mam0QalLXoaLcrcCbJa0jaVVgb5ZMTPVp0myz95OOPIZPJbminJnVWt1PMVUZpJ4LfF3S8cDFEXFtnsejtaLcu0bJfhNpFtb7gB+TKspd1GbfMd99RMzPrz8deBK4BRjuYD4KHB4R50p6D2nq8LexdEW5tYBrJf00zyTbmt06hW7Nh4zMbJDU/TLXxlSUi4jvRcQOEfFm0ljF8D30h7Cks/oflkwj7opyZlZry1IwqB8aUVEut2O4ctyGpCOWM/Om3wK7tbzGcMfhinJmVmuLh6ov/dCIinIR8QRwrqR1chs+3lJb9SPAN/MRydPAcBk+V5Qzs1rr19VJVSlqfg5sOfOHYTZASk7Wt8rRl7Qr87lMXn3q/pW/cxZOPb/r11tWvpPazKxPhup9AOEOwsysX8JzMZmZWTvuIMzMrK2hmg9Su4Mws4G1ytGX9LsJY/IRhJmZteUOwszM2nIHYWZmbbmDMDOzttxBmJlZW0OL6n0V0yCUHN0uP2eupJ9IWiOvd8lRM6u1GIrKSz8MQsnR04Ajc2nR84HP5vUuOWpmtRYRlZd+GISSo5sB1+TH04F358cuOWpmtTYIRxB1Lzl6K7BPfnxAy3qXHDWzWqt7BzEIJUc/BHwrV7q7iHSkAC45amY1V/ermBpfcjQibo+IPSNiR1KVubvzU1xy1MxqbWjRUOWlHxpfcrRl/QTgi8C0/BSXHDWzWqv7KaYqYxDbANdLmg0cBRwz2o6Spkg6DSAiFpPGJmZImksq+9lacvQISXeRxiSeLzkqaSFwBPBFSQuHL1sllRy9DfgJS5ccPUjSAtKX/29JZUYhlRxdnTRGMQuXHDWzmql7B+GSo0vzh2FmVXVdAnSFw3et/J3z3L//3CVHzcxeLOo+SO0OwsysT+o+1YY7CDOzPqn7EcQy3ertZanb3qc2LbtpuU1ssz8LfxaDtJSei2mQTW1gdtNyS2Y3LbdkdtNyS2aXbHPjuIMwM7O23EGYmVlb7iA6V3Jmv1LZTcstmd203JLZTcstme0ZO1v4RjkzM2vLRxBmZtaWOwgzM2vLHYSZmbXlDsLMzNryVBt9lmt070wquxqkKcuvjx5fPSBpdVIZ2Hsi4rEuclYEnhtun6S3kGqE3BYRl9Utd8RrTCGVpF0E3BkRXdUHKd3mXre3JXdD4ImIeEzSZFIhrdsj4tYXU27p7IHQ71u5m7CQ/ic9C7gW+CdghZZtF3SRuydwF3AZcFpeLs/r9uyyzd9pebwrqYDSVcD9wN5d5N4CrJUffxb4BalQ03Tg2Lrl5rzdgBuAnwKPAheTKgz+DNigbm0u1d6cfSRwL6l+yqH5z+8B84AjXiy5pbMHZel7A5qw5P/hDwO2B/4zfxGsk7fd3EXufGBym/UbAfO7bPNNLY+vAnbIjzcGbugi99aWxzcAq+THk4A5dcsd/jsC1m35bM/Pj98OXFm3Npdqb86YB6xCKtT1x5bXWa31/Qx6bunsQVk8BlHNuhExLSJmR8Qnge8A10h6Ld0VGZpEqt090m+AFbrIHWmNiLgJICLuASZ2kfWEpK3z44eAlfPjSXQ3plUqF2BiRDyYH98HvAYgIqaTTu11qlSbS7UXYHFE/Bl4DPgz8HDO/tOLLLd09kDwGEQ1K0haOSKeBoiIH0j6PXAF6beNTp0OzJJ0FunUD6TTWe8ll2HtwuaS5pCqXk2WtFZEPJprd3fT+RwG/FDSLcADwA2Srga2Bf61hrnkrO8BM4B9SadqkLQq3XWWpdpcqr0AN0n6Eenf7QzgDEmXk2vGv4hyS2cPBN9JXYGkw0mnbK4esf71wAkR8fYusrcE9iH9ZijSEcVFEdHVP1BJrxmx6rcR8ZyklwFvjojzuvtgz5kAAArpSURBVMieSBo/2ZQlR0FXRBeD34VzVwA+AmxJGjc4PSIWS1oFWC8ifl2nNhdu7yTgANKR7zmkCyQOJh2pnNTpb89tcv8COKiuuaWzB4U7CLMekrRORDzclFyzsXgMokOS/rcHGS+VdJyk2yU9nJf5ed2avWjnKK/bk8tGe5kr6RWSTpZ0kqR1JH1F0hxJZ0tav8t2DX/O83v5Oefnvyw/niLpHuA6Sb+WtFvdcnPeGpKOlfR/JR08Ytt3ushdXdK/SJon6XFJD0qaKemQLtu7k6SrJP1A0gaSpkt6TNKsfATfTXaRNg8SH0FUkM/lL7WKdErhDoCI2LbD3CuA/wXOiIjf53WvAD4A7NHlqasdRtsEXBwRHX3pFsy9HLiEdD74YOCHwJmkc/Bvi4h9O8nN2aN9zofk7I4+Z0lzI2Kb/Pgq4HMRMUvSpsCPImJKnXJz3rnAncBM4EPAc8DBEfGMpJsiYrS/3/FyLwTOJ12a+x7S3+NZpMt+fxMR/9Rh7vXAPwNrAicAh0fEOZL2AI6JiDd0kluyzYPEHUQFki4CngCOIV3tINI9EbsCdHpOWNIdEbHZsm6rmL0YuDq3daRdImKVmuXeHBGvz4/vi4gNW7bNjojtO8nNzy/yOUu6Hdg6IhZJmhkRu7Rse/5Lvi65+flLfZaSjgL2Jo2DTe+ig7glIrZr+XlWROykdFHEbRGxeYe5Y/27eH5bndo8SHwVUwURsY+k/UlzxX89Ii6S9Fw3g4XZryV9jvSb7R8AJL2cdARx/1hPrGA+8PcRcefIDZK6yS6V23q687/H2NaJUp/zScClko4DLpf0H8B5wB7A7BrmAqwkaUJEDAFExNckLQSuAVbvIvdPknaNiJ9L+mvgkZw/JKndLxNVPS1pT+ClQEjaLyIuyKfaFneRW7LNg2N53GwxKAvpEPTfgIuAhT3IWws4nnQH56N5mZ/Xrd1l9t8Am42ybb8a5v4LsHqb9a8Dzqnx57w78GPSzW1zgUtJdY1XqGnuCaTTaiPX70WazqPT3G2B64HHgZ8P/xsB1gU+1UXudqTLyS8DNge+SbpvYR7wxi4/iyJtHqTFp5g6IGk74A0RMa3fbTEzK8UdRJckbR5dTKImaXPSQGzrZH0XRcT8XrStRHbTcktnj/J6H4yI79cxV9LOQEQa+N6SdPRwe0RcWtPc1wL70zJxIWmw/olucsfIPjMiHu82exD4MtfuXdnpEyV9nnTVhEiHurPy4zMlHdlNo0plNy23dPYYvlrHXEn/DHwLOFnSscC3SWMPR+YB67rlfgqYRprGZCfS3EkbADMl7d5p7jjZv+w2e1D4CKICSd8abRNwSESs0WHuAmCriHhuxPoVgXkRsUknuSWzm5ZbuM0jL39+fhOwaUSsVKfcnD2XNOnkSsDvgVdHxBNKd2lfF51fsl00N9Kd5KsCl0bE7krTdF8Y3V3FVCx7UPgqpmo+CPwj8EybbQd1kTsEvBIYeTXU+nlbN0plNy23ZPbLgb8kDXq3EmnG37rlAiyKiMXAU5LuHj5NExF/ltTNZ1EqF9L31GJS5/OSnHuf0pQk3SqZ3XjuIKqZRZr+9wX/c0r6She5nwZmSLqTJZdbbki6cucTXeSWzG5absnsi0lXXr3g0lNJP6thLsCzklaNiKeAHVtyX0p3nWWp3NNIE1rOBN5MuvIMSeuSL0utafZA8CmmCiStDTyd//H3OnsCSyrKDU/WNyv/NlbL7Kblls5uEkkrRcQLjoSVpvZYPyLm1ik3Z2wFbEH6Ja0nVfWWR/YgcAdRA/mmreevrol8M1eds5uWO8brrR4RT77Yc0tmNy23dHaTuIOoQNIawBeAVwOXRcSPWrZ9JyI+1mHu9qSrKF5K+o1W+TUeAz4WuchPnbKbllvhdZeavuHFmlsyu2m5pbObxGMQ1XyfdH30ucCHJL2bPMEZsMuYzxzbf5GmrbiudaWkXfJrbtfuSX3Oblouko4YbRNdTC/RtNyS2U3LLZ09KHwfRDWvjYgjI+KCiNgHuAn4X0nrdJm72sgvRICImEl3lepKZjctF1J1t7VIV6m0LqvT3f8DTcstmd203NLZA8FHENWUmuDsMkmXkCanay05+n7g8m4aXDC7abmQOvQLIuLGkRskHfoiyi2Z3bTc0tkDwWMQFUg6AbgyIn46Yv1ewH92eRPXO1gyBURrydGupicomd3A3M2ARyLiwTbbXt7pQHjTcktmNy23dPagcAdhZmZt+RRTRSow2ZtS0fQPA/uNyL0Q+N7IqSHqkN203BHZ+5PuqO51mxuR28Q2N/GzGCQ+gqhAabK3g0gTvi3Mq18NvBc4KyKO6zD3TNJlnGeMyD2EVKfgwC7aXCS7ablNbLM/i/K5pbMHhTuIClRusrexSmEuiIhNO8ktmd203JLZTcstmd203NLZg8KXclUzPNnbSN1O9vaopAPyNBAASJog6UBeOFFbXbKbllsyu2m5JbOblls6ezBEDcra1X0hFT65i1T28NS8XJ7X7dVF7mRSWckHgQWkm/EezOs26rLNRbKbltvENvuzaPZnMUiLTzFVpMKTveWb7hQRD/Uib3lkNy23ZHbTcktmNy23dHaTuYNYBurxRHJKhUkeiIinJQn4ALADcBvw3YhYVLfspuU2sc3+LMrnls4eFB6DqEDS9kpzxv+MNGf8icDVkmZK2qGL6EtZ8ndwHPBO4DpS+cNTu8gtmd203JLZTcstmd203NLZg6Hf57iasACzgb9os34X4JYucm9reXwjMKHl545zS2Y3LbeJbfZn0ezPYpAWH0FUU2oiufslvTU//hVp7qHh86HdKpXdtNyS2U3LLZndtNzS2QPBYxAVSPoW8FraTyR3b0R0VLZS0gY5cyLwOLArcDNphsnPRMSMLtpcJLtpuU1ssz+L8rmlsweFO4iKVHZSvS2ATUlTnwxfHdVtofei2U3LLZndtNyS2U3LLZ3ddO4gzMysLY9BdEnS1EK5xa6iKJXdtNyS2U3LLZndtNzS2U3iDqJ7KpR7SqHcktlNyy2Z3bTcktlNyy2d3Rg+xWRmZm25HsQykLSYdJPcFyL3rJJuiohubpZD0qbAZ4HX0PJ3EhFvHfVJfc5uWm7J7KbllsxuWm7p7KbzEcQykDSHNEnf64EDI+IRSTdHxOu7zL0FmEa6Wef5uZ2iTa3cumQ3LbdkdtNyS2Y3Lbd0dtP5CGLZLIqIz0l6D3CtpPeT5mXqRe7JPchZntlNyy2Z3bTcktlNyy2d3Wg+glgGrUcLkrYkVZjbMCLW7DBv7fzwU8ADwPnAM8PbI+KRLtpaJLtpuSWzm5ZbMrtpuaWzB4U7iGUg6XPAyRHxR0lfBN4KXBURR3eYdy/pCKT1Sqjn/0IiYuMu2loku2m5JbObllsyu2m5pbMHxvKe/KnJCzAn/7krcC3pzurrepD7HmCN/PhLpN9kduhRm4tkNy23iW32Z9Hsz2IQlr43oEkLcHP+81jg4Pz4ph7ktnY81/Sq4ymZ3bTcJrbZn0WzP4tBWHyj3LL5jaRTSL9xXCppJdJEX90avnLincC0iLgQWLEHuSWzm5ZbMrtpuSWzm5ZbOrvR3EEsm/cAV5DqUD8GrE26frpb7TqeXv3dlMpuWm7J7KbllsxuWm7p7Gbr9yGMlwBYFXgXsEn+eX1gzzpnNy23iW32Z9Hsz2IQFl/FZGZmbfkwyszM2nIHYWZmbbmDMDOzttxBmJlZW/8fl2Iy6u+bUg8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "df = pd.DataFrame(pairs, index=range(0,len(pairs)), columns=list(['Name1','Name2','pvalue']))\n",
    "#pvalue越小表示相关性越大,按pvalue升序排名就是获取相关性从大到小的股票对\n",
    "df.sort_values(by='pvalue')"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Name1</th>\n",
       "      <th>Name2</th>\n",
       "      <th>pvalue</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>sh.601328</td>\n",
       "      <td>sh.601988</td>\n",
       "      <td>0.015970</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>sh.601328</td>\n",
       "      <td>sh.601998</td>\n",
       "      <td>0.021157</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>sh.601009</td>\n",
       "      <td>sh.601169</td>\n",
       "      <td>0.023653</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>sh.601988</td>\n",
       "      <td>sh.601998</td>\n",
       "      <td>0.047497</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Name1      Name2    pvalue\n",
       "1  sh.601328  sh.601988  0.015970\n",
       "2  sh.601328  sh.601998  0.021157\n",
       "0  sh.601009  sh.601169  0.023653\n",
       "3  sh.601988  sh.601998  0.047497"
      ]
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "source": [
    "# plot(prices_df[['sh.601328','sh.601998']], chart_type='line', title='Price')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "source": [
    "# ols\r\n",
    "x = prices_df['sh.601328']\r\n",
    "y = prices_df['sh.601998']\r\n",
    "X = sm.add_constant(x)\r\n",
    "result = (sm.OLS(y,X)).fit()\r\n",
    "print(result.summary())"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              sh.601998   R-squared:                       0.682\n",
      "Model:                            OLS   Adj. R-squared:                  0.682\n",
      "Method:                 Least Squares   F-statistic:                     1323.\n",
      "Date:                Thu, 09 Sep 2021   Prob (F-statistic):          1.20e-155\n",
      "Time:                        16:01:31   Log-Likelihood:                -566.43\n",
      "No. Observations:                 619   AIC:                             1137.\n",
      "Df Residuals:                     617   BIC:                             1146.\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.3818      0.226      1.687      0.092      -0.063       0.826\n",
      "sh.601328      0.8602      0.024     36.378      0.000       0.814       0.907\n",
      "==============================================================================\n",
      "Omnibus:                        0.497   Durbin-Watson:                   0.070\n",
      "Prob(Omnibus):                  0.780   Jarque-Bera (JB):                0.340\n",
      "Skew:                           0.003   Prob(JB):                        0.844\n",
      "Kurtosis:                       3.115   Cond. No.                         90.0\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "source": [
    "import matplotlib.pyplot as plt\r\n",
    "fig, ax = plt.subplots(figsize=(8,6))\r\n",
    "ax.plot(x, y, 'o', label=\"data\")\r\n",
    "ax.plot(x, result.fittedvalues, 'r', label=\"OLS\")\r\n",
    "ax.legend(loc='best')"
   ],
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1982554c640>"
      ]
     },
     "metadata": {},
     "execution_count": 8
    },
    {
     "output_type": "display_data",
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFlCAYAAAAzqTv+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3xT5f0H8M/TECBFsQXRQbiOaVFEWqnKxBui4oVLhSkydF42UeemoFZhOIF5AayKTn9u4mW4CQoIVAEVnOBQJs5iW6EKXhEIKhUogg2Qts/vj9OTpMk5yUlyTk4un/fr5cueJE2eIOab7/N8n+8jpJQgIiKi5MqxewBERETZiAGYiIjIBgzARERENmAAJiIisgEDMBERkQ0YgImIiGzQKpkvdvTRR8uePXsm8yWJiIhss2HDhh+klJ207ktqAO7ZsycqKiqS+ZJERES2EUJ8o3cfp6CJiIhswABMRERkAwZgIiIiGyR1DViLz+fDjh07cPDgQbuHYpm2bduia9eucDqddg+FiIhShO0BeMeOHTjyyCPRs2dPCCHsHo7ppJTYvXs3duzYgV69etk9HCIiShG2T0EfPHgQHTt2zMjgCwBCCHTs2DGjM3wiIoqd7QEYQMYGX1Wmvz8iIopdSgTgVDJt2jQ8/PDDuveXl5fjk08+SeKIiIgoE6VdAC6v9GDQzNXoNWkFBs1cjfJKT3JfnwGYiIhMkFYBuLzSg8lLNsJT54UE4KnzYvKSjQkH4QceeAAFBQU4//zzsWXLFgDAM888g1NPPRX9+/fH6NGjUV9fj//+97947bXXUFpaisLCQnz55ZeajyOizGV3EkCZI60CcNnKLfD6Glvc5vU1omzllrifc8OGDXj55ZdRWVmJJUuW4MMPPwQAjBo1Ch9++CGqq6txwgkn4LnnnsMZZ5yBESNGoKysDFVVVejdu7fm44goM1mVBFB2ihqAhRDPCyF2CSE2Bd12nxDiYyFElRBilRCii7XDVOys88Z0uxHvvvsuLrvsMuTm5qJ9+/YYMWIEAGDTpk0466yz0K9fP8ybNw81NTWav2/0cUSU/qxIAih7GcmA5wK4KOS2MinlyVLKQgDLAdxr9sC0dMlzxXS7UVpVytdeey2efPJJbNy4EVOnTtXdRmT0cUSU/qxIAih7RQ3AUsq1APaE3PZj0GU7ANLkcWkqHVoAl9PR4jaX04HSoQVxP+fZZ5+NpUuXwuv1Yv/+/Vi2bBkAYP/+/ejcuTN8Ph/mzZvnf/yRRx6J/fv3+6/1HkdEmceqJICyU9xrwEKIB4QQ2wGMQ4QMWAgxXghRIYSoqK2tjfflAAAlRW7MGNUP7jwXBAB3ngszRvVDSZE77uc85ZRTMGbMGBQWFmL06NE466yzAAD33XcfTj/9dFxwwQXo06eP//FXXnklysrKUFRUhC+//FL3cUSUeaxIAih7CSmjJ69CiJ4AlkspT9K4bzKAtlLKqdGep7i4WIaeB/zpp5/ihBNOMDretJUt75Mo05VXelC2cgt21nnRJc+F0qEFCSUBlNmEEBuklMVa95nRC3o+gBUAogZgIqJ0V1LkZsAlU8Q1BS2EOC7ocgSAzeYMh4iIKDtEzYCFEC8BOBfA0UKIHVAy3UuEEAUAmgB8A+AmKwdJRESUaaIGYCnlWI2b2W2CiIgoAWnVCYuIiChTMAATERHZgAEYwI4dOzBy5Egcd9xx6N27N2677TYcPnwY77zzDoYNGxb2+OXLl6OoqAj9+/fHiSeeiKefftqGURMRUTrL+gAspcSoUaNQUlKCzz//HJ999hkOHDiAKVOmaD7e5/Nh/PjxWLZsGaqrq1FZWYlzzz03uYMmIqK0Z8Y+4LS2evVqtG3bFtdddx0AwOFwYPbs2ejVqxcGDx4c9vj9+/ejoaEBHTt2BAC0adMGBQXsgkNERLFJrQA8YQJQVWXucxYWAo89pnt3TU0NBgwY0OK29u3bo3v37vjiiy/CHt+hQweMGDECPXr0wJAhQzBs2DCMHTsWOTlZP5lAREQxyPqoIaXUPA1J73YAePbZZ/H222/jtNNOw8MPP4zrr7/e6mESEVGGSa0MOEKmapW+ffti8eLFLW778ccfsX37dvTu3Vv39/r164d+/frh6quvRq9evTB37lyLR0pERJY5fBi49FLguuuAX/86KS+Z9RnwkCFDUF9fj3/+858AgMbGRtxxxx249tprkZubG/b4AwcO4J133vFfV1VVoUePHskaLhEZUF7pwaCZq9Fr0goMmrka5ZUeUx5LGWrZMqBNG+Df/waefDJpL5v1AVgIgaVLl2LRokU47rjjcPzxx6Nt27Z48MEHAQBvv/02unbt6v+nsrISDz30EAoKClBYWIipU6cy+yVKIeWVHkxeshGeOi8kAE+dF5OXbNQMrLE8ljKQzwd07w6MGKFcX3opsG5d0l4+taagbdKtWzcsW7Ys7PZzzz0XXq837Hb1zGAiSj1lK7fA62tscZvX14iylVvCTjGK5bGUYd54A7jkksD1Rx8BRUVJHQIDMBFllJ114V+a9W6P5bGUIRoagIIC4KuvlOsLLgBWrgR0im6tlPVT0ESUWbrkuQzfHstjKQO89RbgdAaC7//+B6xaZUvwBRiAiSjDlA4tgMvpaHGby+lA6dDwhjmxPJbSWGMjcOKJwIUXKtfnnAM0NQGnnmrrsFJiCjrSnttMIKW0ewhEWUNduy1buQU767zokudC6dCCsDXd8kqPfw3YIQQapYRb57GUxlavBoYMCVz/97/AL39p33iC2B6A27Zti927d6Njx44ZGYSllNi9ezfatm1r91CIskZJkTtiEFWrn9UCrEYp/Zkvg2+GaGoCTjkFqK5WrgcOVCqcU6hroe0BuGvXrtixYwdqa2vtHopl2rZti65du9o9DCJqxurnDPfuu8DZZ7e8PvNM+8ajw/YA7HQ60atXL7uHQURZhNXPGaqpScl0P/xQuT7lFOXnFMp6g6XmqIiILMTq5wz0/vuAwxEIvmvWABs2pGzwBRiAiSgLaVU/CwCD+3SyZ0BgS8y4SQmcdRZwxhnKdd++yl7fNDinnQGYiLJOSZEbowe4EVz2KQEs3uCxJfDF0xKTARvKPt6cHOC995Trt94CNm1SMuE0wABMRFlpzeZahG4QVAuxki1SUZiWrO9hLSVw3nnA6acr18cdp/R1Pv98e8cVIwZgIspKqVSIFetYYg3YGUVd112zRrl+803gs8+AVrbXFMeMAZiIslIqFWLFOpZU+vKQNFICF18MFBcr1z16KGf4Dh1q77gSwABMRFkpldpQxjqWVPrykBTV1UrW++abyvWyZcDWrUpf5zTGAExEWamkyI0Zo/rBneeCAODOc2HGqH62NOKIdSyp9OXBUlICI0cChYXKdefOStY7bJi94zKJSGaf4uLiYllRUZG01yMiylRqL+tI/a7T2qZNQL9+geulS4GSEvvGEychxAYpZbHWfem3ak1ERFH7Xae1K64AFi1Sfu7YEfB4gDZt7B2TBTgFTUREqeHTT5WzedXgu3Ah8MMPGRl8AWbARESUCq66Cpg3T/n5iCOA2logw0+RYwZMRET2+ewzJetVg+/8+cD+/RkffAFmwEREZJfrrwf+8Q/lZ6cT2LcPcGXoVioNzICJiCi5vvxSyXrV4PvCC8r2oiwKvgAzYCIiSqabbwb+/vfA9YEDQLt29o3HRsyAiYjIelu3KlmvGnyffVZptJGlwRdgBkxEaSDjm05kultvBZ54InC9f79S6ZzlGICJKKWpR++pp/+oR+8BYBBOddu2KYcmqP7+d+DGG+0bT4rhFDQRpbSsPnovnd15Z8vgu28fg28IBmAiSmlZefReOvN4lLXeRx5Rrp94Qlnrbd/e3nGloKgBWAjxvBBilxBiU9BtZUKIzUKIj4UQS4UQedYOk4iyVdYdvZfOJk8GunYNXO/dC/zhD/aNJ8UZyYDnArgo5La3AJwkpTwZwGcAJps8LiIiAFl09F46++47JeudOVO5fvRRJevNY24WSdQALKVcC2BPyG2rpJQNzZfrAXQN+0UiIhOk0rm9pGHqVOWcXtXu3cDEifaNJ42YsQZ8PYA39O4UQowXQlQIISpqa2tNeDkiyjYlRW6sm3QeZo9RDmafuKAKg2auRnmlx+aRZbFdu5Ss9y9/Ua5nzVKy3g4d7B1XGkloG5IQYgqABgDz9B4jpZwDYA4AFBcXy0Rej4iyl1XbkbjHOA733w/8+c+B69pa4Oij7RtPmoo7AxZCXANgGIBxUkoGViKylBXbkdSg7qnzQiIQ1JlZ6/jhByXrVYPv/fcrWS+Db1ziyoCFEBcBuBvAOVLKenOHRETZKFomasV2pEhBnVlwiIceAu6+O3D9/ffAMcfYN54MYGQb0ksA3gdQIITYIYT4LYAnARwJ4C0hRJUQ4u8Rn4SIKAIjmagV25G4x9iAPXuUrFcNvvfeq2S9DL4JM1IFPVZK2VlK6ZRSdpVSPiel/IWUspuUsrD5n5uSMVgiykxGppet2I7EPcZRzJ4NdOwYuP72W2D6dPvGk2HYCYuIbGckE7ViOxL3GOuoq1Oy3ttvV64nTVKy3p/9zN5xZRgexkBEtuuS54JHIwiHZqIlRW5T12bV52IVdJAnnwT++MfA9Y4dgDuL/zwsxABMRLYrHVrQYosRkLxM1OygnrZ+/BE46qjA9e23B/o5kyUYgInIdsxEbfb008BNQaU833wDdO9u33iyBAMwEaUEZqI2OHAAOPLIwPUf/wj89a/2jSfLsAiLiCgbPfdcy+D79dcMvknGDJiIko7tH23000/AEUcErm+8Efg7WznYgQGYiJLKqp7Odki7LxL/+hfwm98Err/4Aujd277xZDlOQRNRUlnR09kOadVH2usFWrcOBN/rrlP29TL42ooBmIiSKlPaP6bNF4mXXwZycwGfT7nesgV4/nl7x0QAGICJKMkypf1jyn+ROHhQWesdO1a5HjdOyXqPP97ecZEfAzARJVWmtH9M6S8Sr7wCuFxKwRUAfPIJ8OKL9o6JwjAAE1FSWdHT2Q4p+UXi0CGgQwfg8suV68svV7LeE06wb0yki1XQRJR0mdB0I+W6d5WXA5ddFrjeuBE46aTkDiHdqsJtxgBMRBSnlPgicfgw0KMH8N13yvXIkcDSpcppRkmUSdvLkoVT0ERE6Wr5cqBNm0DwrapSMuEkB18gjarCUwgzYCKidOPzAb/4BbBtm3J90UXA66/bEnhVKV8VnoKYARMRpZM331SaaqjBt6ICeOMNW4MvkOJV4SmKAZiIKB00NCh7eC++WLkeMgRoagIGDLB3XM1Ssio8xXEKmogo1f3738AFFwSuP/gAOO00+8ajIeWqwtMAAzARUapqbAT69wdqapTrM88E1q61fbpZT0pUhacRTkETEaWi//wHaNUqEHzXrQPefTdlgy/FjhkwEVEqaWoCiouBykrl+rTTgPffB3KYL2Ua/hclIkoV770HOByB4Lt2rbLey+CbkZgBExHZrakJGDQIWL9eue7fH/joIwbeDMf/ukREdlq/Xsl61eC7erXS0YrBN+MxAyYisoOUwDnnKIVVgHJi0caNSjCmrMCvWEREyfbhh0qGqwbfVauUM3sZfLMKM2AiomSRErjwQqWxBgD07g1s3qxsN6KswwyYiCgZKiuVrFcNvq+/DnzxBYNvFuN/eSIiK0kJDBumBFwA6NoV+OorwOm0d1xkO2bARERW+fhjJetVg++rrwLbtzP4EgBmwERE5pMSGD0aWLpUuT72WOX4wNat7R0XpRRmwEREZqqpUbJeNfguXgx89x2DL4VhBkxEZJaxY4GXX1Z+zstTAm+bNvaOiVIWM2AiokRt3qycUqQG3wULgL17GXwpImbARESJuOYa4J//VH7OzQV27wbatrV3TJQWmAETEcXj88+VrFcNvi++CPz0E4MvGRY1AAshnhdC7BJCbAq67XIhRI0QokkIUWztEImIUswNNwDHH6/87HAogXfcOHvHRGnHSAY8F8BFIbdtAjAKwFqzB0RElLK++krJep99Vrn+xz+AhgZl6pkoRlHXgKWUa4UQPUNu+xQAhBDWjIqIKNXccgvw1FOB6wMHgHbt7BsPpT2uARMRRfLNN0rWqwbfOXOURhsMvpQgywOwEGK8EKJCCFFRW1tr9csREZlnwgSgZ8/A9Y8/Kuu/RCawPABLKedIKYullMWdOnWy+uWIiBK3fbuS9T7+uHL91FNK1nvkkfaOizIK9wETEQW76y6grCxwXVcHHHWUfeOhjGVkG9JLAN4HUCCE2CGE+K0Q4jIhxA4AvwSwQgix0uqBEhFZaudOJetVg+/jjytZL4MvWcRIFfRYnbuWmjwWIiJ7TJkCPPhg4HrPHiA/377xUFbgFDQRZa/vvgM6dw5cP/wwcMcd9o0nC5RXelC2cgt21nnRJc+F0qEFKCly2z0sWzAAE1F2mj4dmDYtcP3DD0DHjrYNJxuUV3oweclGeH2NAABPnReTl2wEgKwMwtwHTETZpbZWWetVg++MGcpaL4Ov5cpWbvEHX5XX14iylVtsGpG9mAETUfZ48EFlvVe1axfA7ZGWCZ1u9tR5NR+3U+f2ZLJjapwBmIgy3+7dwNFHB66nTwfuvde+8WQBrelmAUBqPLZLniupYwtl19Q4AzARZbaHHwZKSwPX330HHHusfePJUKEZ5E+HGsKmmyUQFoRdTgdKhxYkc6hhIk2NMwATEcVq716gQ4fA9ZQpwP332zeeDKaVQeqRANx5rpSqgtabArd6apwBmIgyz+OPK32cVTt3ttxuRKbSyiD1uPNcWDfpPItHFBu99Wmrp8ZZBU1EmWPfPqXCWQ2+d92lVDgz+FrKaKaYCtPNWkqHFsDldLS4LRljZQAmoszw1FNAXl7gevt2YNYs+8aTRfQyxfxcJ9x5Lggome+MUf1sn27WUlLkxoxR/ZI+Vk5BE1F6278faN8+cD1xIvDoo/aNJwuVDi1osQYMKMVWe+t9yG3dCrPHFKZk4A1WUuRO+hgZgIkofT3zDDB+fOB661agRw/bhpOt1MBVtnJL2HajbO92FQmnoIko/Rw4oKz1qsH3lluUtV4GX9uUFLmxbtJ5cOe5wvb6ZnO3q0gYgIkovcydCxx5ZOD6q6+AJ5+0bTjUkl1betIRp6CJLMJTX0xWXw8ccYSS6QLADTcAc+bYOyYKk4wtPZny/xYzYCILqI0JPHVeSCjrYBMXVKHnpBUYNHM1yis9dg8xvcybB7RrFwi+n3/O4JuirN7So/X/1uQlG9Py/ylmwEQW0GpMoFWUoj423b/JW8brBfLzgUOHlOvf/AZ44QV7x5RhzM4mgwuyrPh7bVfbSCswABNZINp6l9fXiGmv1eBQQxPPRtWzYAFw5ZWB682bgYLUa+KQzqw6hMDKLT2ZtMbMKWgiCxhZ76rz+ng2qpZDh5R9vWrwHTtWmXpm8DVdOp7Pq/f/lt0nKsWDAZjIAqVDCyDi/N10/CZvmsWLgbZtleYaAFBTA8yfb++YMlg6ZpN2tY20AgMwkQVKitya556qXE4H8nOdmvdJIPsKtQ4fBjp1An71K+V69Ggl6z3xRHvHleHSMZu0q22kFbgGTGQytahFj0MIzBjVDwDC2vepsmo9+LXXgJEjA9cffwz062ffeLKIVgvJdMgm7WgbaQUGYCIThRa1hHI5HWHf1tX2faHStbLTMJ8P6NlTOSoQAIYNU4KxiHfynmJlRsVypuzJtYOQMtJEmbmKi4tlRUVF0l6PyGqhHz71hxuwt96n+Vh3hA+nXpNWaE5ZCwBfz7zU3EGngtdfBy4Nel8ffQQUFdk3HoqL1hdOrS+Z2UwIsUFKWax1HzNgojhpbeHQI4CIh5DbdSB40jU0AMcfD3z9tXI9dCjwxhvMetNUJu3JtQMDMFGQWKbTtD589EQLpOm6FheTVauUgKv68EOgWDMxoAhSaco3HauoUwkDMFEzo00J1A/ASBlvMCOBtKTIjYpv9uClD7ajUUo4hMDoAZlRaILGRqWa+bPPlOtzzgHWrEnLrNfu4GdV44x4Zc3MjUW4DYmomZGmBMF9aPXkuZwxb5Eor/Rg8QYPGptrMhqlxOINnvTfirR6NdCqVSD4vv8+8M47aRt87e5BnGqNMzJpT64dmAETNTMynRZt2lkAmDaib8zZSMatpTU2KkVVG5t7Xp9xBvDee2kZeFWp8N8oWVO+RjN9q/s+ZzoGYKJmRqbTok07S8Q3FZhRa2lr1yrTzKr33gMGDbJvPCZJhf9GyTrqL5Zp7kzZk2sHTkETNYs2nVZe6YnaXtId5wdhOnYkCtPUBJx6aiD4DhigZMIZEHyB1PhvlIwp31Sb5s5kDMBEzSK1uCuv9OCOhdVR20vG+0GY9mtp69YBDgeg7vN/5x3l55zM+YhJhf9GyWjDmAqZfrbgFDRRkNA1rbKVW1DxzZ4WBVJaIjXZiOd102YtTUolw33/feW6Xz+gslIJxhlGrVSf98E2qH8VRMSvZNaNw8q/F6xsTh4GYKIgWutf89Zvi/gx685zRWyyYVTaraV98AEwcGDg+t//BoYMsW88Fiuv9GDBh9sR/D2s3teE0kXVADKnZ3dW7ElPEQzAREG01r+iTTsP7tMJg2auTq/MNRFSAuedp0wzA8o5vTU1GZn1BitbuQW+xvC/Db4mmb7V6hrSdjYmDTEAEwWJdZ1LQGLBh9v9H8x2N0awXEWFUmilevPNlt2tMlikvxuZtj6adrMxaYoBmCiI3vqXgHYmXO9rCrstrffv6pESuOgipZ0koJxi9PnnSpONDBa8HzZHCN06AK6PUjwy+/8eIh1ajQYqvtmjGXxdTgdGD3BjzeZaw+0nPXVeDJq5OjOm7qqqWp5UtHx5y5OMMlRoPYBe8HXmCK6PUlwYgCkjRerko1VodfvCKjTpLPaOHuDG/SX9cE/5Rry4fpvhMaT9dLSUwMiRwLJlyrXbrZxi5HTaO64k0et6liPg/7uS53LG1fmMCDAQgIUQzwMYBmCXlPKk5ts6AFgAoCeArQCukFLutW6YRMZpBdjSRdWYvqwGdfU+zalEveALAGs216K80oN5MQRfVdpOR2/cCJx8sv/y7qv/gl9OvB4lWRJ8Af11XSmBrZl4RjMlnZFd8nMBXBRy2yQAb0spjwPwdvM1UUrQylx8TRJ7632Q0J9K1KPuB453x6fWB3l5pQeDZq5Gr0krMGjm6hYN/SPdlxS/+pU/+O5xtcdxdy7Fgi6nJP3gAbulQucrymxRM2Ap5VohRM+Qm0cCOLf55xcAvAPgbhPHRRQ3sytSu+S5EnrO0A/sSL12Adh33NwnnwB9+/ovbyqZjDcLAm0k0zabjxP3w5LV4l0DPlZK+S0ASCm/FUIco/dAIcR4AOMBoHv37nG+HFE4vXVevUrmeDiaC2xiOf83mNYHdrReu7acuDNuHDB/vvJz+/YouGEuDrVqHfawTNtuEwn3w5LVLC/CklLOATAHAIqLi5Pft40yUqQsUitziSbXmQMhBH46HPiddq0deOCyQJ/dWJ5TAJof2OWVHt1Abss+088+UxppqF56CbjyShw9czXbEYL7Ycla8Qbg74UQnZuz384Adpk5KKJoImWRalvIOxZWG17vDd3P63K2DL7qv6cvq8Heel/E59JrTal+adCjBrekBb7rrgPmzlV+btMG2LsXcCmvw+lXIuvFe1TJawCuaf75GgCvmjMcImOindhSUuRGU4zFVsG0jl8rKXIjt3Xk76x6QUo9TUkvgxZQAu9PhxrgdLQ89ND0wPfFF4AQgeD7r38BBw/6gy+QnFN3iLKdkW1IL0EpuDpaCLEDwFQAMwEsFEL8FsA2AJdbOUiiUNFObCmv9Oh2LtLrahVKK8hHmwrWClJq5hspG1fvqfP64MwRyM91oq7eZ/664403AnPmKD8LARw4AOTmaj6U069E1jJSBT1W567MPfaEUl6kKdJIAS+4q5VaWFN/uEFzWllr2jdSgZc7z6UZsPQaOujxNUnktm6FynsvNPw7UX39NfDznweun39emYImItuwExalpUgVqoNmrtYMeA4hImaowb+jTgmHtpMsHVqA0leqw07FidSOMJ4CqkSKrkKrw/9V9S/8fMHcwAP27weOOCLu5yciczAAU9rSmyLVC15NUmo+PjiYe+q8LaaoQ/fhahVjRWtHqJc1O4TAkW1boc5rLPs2IvjLRJcfd2HdrOsDdz79NDB+fFzPS0TmYwCmjBNtfThYcLZ4lMsJIYDQmevQfbixro3qTZfPGNVPuT8ko3Y64m/ur0533/P2M/hdRaA28qQJC3HUno4orfRk5bpupN7gmfSalF4YgCktxPJhZnQLzT3lGzFv/bYWBVB6EpkSjjRdXl7pCa8IS2C3fNO27dj6t2v913++4Cb865RhAIAD6X44RJwi7Rm36s/Bjtek9BPvNiSipFE/zDx1XkgEPsz0+hIb2UKjHq5gNNZZ1YCibOUW+EJOgvA1ybAtUIZMmoT3g4Lvybe97A++Kq3tVZkuWuexTHlNSj/MgCnlRfow08smok0Tx3K4QqL7cLWyoYkLqlDxzZ6o+5kN+fZboEsX/+WMC8bj6VNG6D48E9tJRpohMeXPOEZ2vCalHwZgSnl6236M9mbW+nCO5YMwOHPRCurRpse1vkBIAPPWb8NRLmdiRVj33gvcd1/gevdunPCNF+4IvaszrZ1ktOneWGoCzGLHa1L6YQCmlOfQaajhEELj0S2FrvOqH856gU+Pp86L0leqMe21GuzzBhpkAMDtC6rQFPS42xdUAUDUDExC6YXhcjpib/n4/ffAz34WuH7oIaC0VHndDvCvL2dDO8loMySD+3QKW26w+s+BrTzJCAZgShl6maReB6lofZ711nm9vka0deaEBb5ofI3SH7Q9dV5MaA60oZoAlC6q8gfgSM076up9mD2mMLZq2fvuUzJfVW0tcPTRYQ/LltN8Ik33lld6sHiDp8XfAQFg9ABru3xly589JUbIBPrlxqq4uFhWVFQk7fUotYVuAfrpcEOL7TjqVh29owAdQqBJyrAPN/V5o01RP9Yc+Mw6ujCUuj8YACYuqNJcc9Y7uJtaI6AAACAASURBVEFTbS1wTNDJnw88APzpT4kPNM0N0jm5yR3hcIuY/tyJEiCE2CClLNa6j1XQZIvQyuY6ry+su5Q6jVg6tAAupyPsORqlDKuKDn7eSNTJ63WTzsNjYwoRfTI7dnVen38tctzA7mGvEdOU5MyZLYPv998z+DbT+vuh/tmyGIpSGQMw2cJof+Sddd6wbUVaa79qsDb6vLJ5DIAyXXhG7w6xvgVD1HHdX9IPs8cUxn660O7dykLx5MnK9dSpSqeQ4GCc5SJtO9MremIxFKUCrgGTLYxmIOoHZfC2ol6TViT0nFqP37rb+i0pMZ8u9OijwB13BK6//bZl4RX56f3ZshiKUhkzYLKFkQxE74MyUlaTl+uMawxGgnd+rhP5MTy/1usYUlenZL1q8P3Tn5Ssl8E3ZjzXmFIZM2CyhVZmEkxA+2zdSL/706EGHG4wXtUcHNwjVSoDypeBqcP7+rf3THutxtA2ppizrSeeAG69NXDt8bRoskGx47nGlKoYgCnp1CrlSGu1kWrztU4kAiL3cg6V53K2qJquP9wQ9hj1VCR3SJW1+oGutW0KiHPryb59QF5e4PrOO4GyMsPvh4jSD7chUVJpNYfQE22riN72k2jU7U16zSqA6EcMmurpp4Gbbgpcb9sGdOtm/esSkeW4DYlShtEqZSDyumx5pcdw8HXmCOTnOjXXAPXG065NqxYZ8qCZq9Fr0goMmrla9xCImB06BEyZEgi+t96qrPUy+BJlBU5BU1LFkrHqFS+VV3pQuqja0HMETx+rU8YTF1ShbOUWDO7TKWKfaTXQWnKs3EcfAddcA2zaBFx7LfDgg0DnzvE/XwrjubhE2hiAyXLBH8BGRSpe0jrCT0vwFLZWw/4X12+L+Pulr1SjXetWMZ/EFNHhw0oHqwceUPbyLl8OXHpp7M+TJnguLpE+BmCyVOhhCNEIIGqWZDSQBwfwWKa+VcG9n+MdQwtVVUq2W10NXH018PjjQH5+7M+TRuI5SpIoWzAAkyXKKz1hVcrR5Oc6UXnvhWHPEzp9GW3LkBazWw/GtLfX5wNmzFAOUTj6aODVV4ER+uf1ZhK2giTSxyIsMp067RhL8AWU+iOt51H7RavTl4P7dDL0fNNeq/H/HEuDjmD5uU7dPsOGbNwInH660kLyiiuUNd8sCb5A5KYp2cSyQj5KawzAZLrpy2pinu4FgH0h071605fR1m5VwdPH8ey2czoEpg7vG18npYYGZZ13wAClmcaSJcC8eUDHjrEPJI1FOighW+h9kWQQJk5Bk6nuKY8981WFZkVmTlOGBvdQOQIIruvKz3X6O18BMRYM1dQoFc4bNgBjxgBPPql5Xq8qk6uEeS4u18FJHwMwmaa80oN5BrPTUFpZUTxrvcGC+zbrPVdopbQaKHJbx/G/RkMD8PDDynRz+/bAokXAr34V8VeyoUo421tBch2c9HAKmkxTtnKLoWpnd54LVw3sHnVat3RoQULn9EoJ/zTf4D6dNJ/LU+dF0V9W4Z7yjYlNE376KTBokHJs4PDhShYcJfgCkbMjygxcByc9zIApLlrTptG+0QsAX8+Mbc9rIo1S67w+lC6qxpSlG/HTYf016b31Ps11Za+vEXcsVBp+6GZwhw8rnat27VLWd19+WSm20jizWAuzo8zHIxFJDwMwxUxv2jQv1xlx/TeWb/zqayTK1yThixB8o2mUUn9KeNmylhXNmzbFfGTgUS6n5l5jZkeZg+vgpIcBmGKmN216MELlc6zf+ONpnGGVsIIZnw/o3RvYvl25vvRSJRgbzHpV5ZUe/KRxCpMzR6RVdpTJRWRmyfZ1cNLGAEwx05seDZ0uViuLQ4/zCxXPdHay+cfzxhvAJZcE7vjoI6CoKK7nLFu5Bb7G8En2I9q2SpsP62woIiOyCgMwGaYGSqPrsu3bOlE19cKIj9H6AJ+woEpJJpN3UmZU3Y5srWS9X32l3HD++cCqVTFnvcH0vmTUxbmNyw7cYkMUP1ZBkyHBzQSMqvP6olYR6001W3lMda4z8l97p6NlUB2yvRpr77kwEHz/9z/grbcSCr5AZlTHsoiMKH4MwGRIvGuy0bbTJPuD+qqB3ZHfrk3kB0llD7GjqRFr/vF7PDd/inL7OecATU3AqaeaMpZM6BKVCV8iiOzCAEyGxBsoPXVe9Jy0Ar0nv457yltWNZdXepCTYBYZizyXE4s3eKJm8b4miTO3b8SXZSPRa1fz9qT//hd4552Es95gJUXu+NpcppBM+BJBZBchrZzrC1FcXCwrKiqS9npknkEzVyfUlSqYO8+Fnh1dWPflHlOeL1S71o6wfb8upwNtnTlR22QK2YQVc2/Dibu+Vm4YOBBYtw7Iyc7vqkYqnFkFTaRPCLFBSlmseR8DMBlRXunBhAVVdg8DgJLJ7jvo01wnVltLlld6MO21Gv8eW62gHOrU7ZuwaP6kwA3vvguceaaZQ08roQVygPJFJt2ydCI7RQrA2fm1nmJWUuRGnkv7SL88lxOPjSlEjkmzswKAQ2eqN8/lxKGGJs3gGzr1eaihyf9zpOArZBPK/znRH3xrfvYLlFdsy+rgC7BNJpHVGIDJsGkj+mqu900b0RfTl9W0OE0oEV3yXGjUmZmp8/o0i8EcQrTIzIwWjZ3i+RRfPzQChd9+DgD4ww2P4PPX30HJgG4JvIPMwApnImsltA9YCHEbgBugJC3PSCkfM2VUlJIitdQzc3o6nrXmJilbTItGDRJSYuH8u3Hajk8AAFuO7o7fTXwW7/7pgphfO1PpnSDFCmcic8SdAQshToISfE8D0B/AMCHEcWYNjKxRXunBoJmr0WvSCgyauTqtDgV3OR0tjhgMFhoUIgWJ/ju3YOtDw/3Bd9yY+zH0t0/hnBNj6+Oc6VjhTGStRDLgEwCsl1LWA4AQ4j8ALgPwkBkDI/Ml2jYw0u/n6RwqYBa1nSUAzcKgwX06YdDM1f7MPLe1xndLKTFvwRQM+uZjAMBX+V1wwe/+hsYcJcis2Vxr2fjTEQ8RILJWIgF4E4AHhBAdAXgBXAIgrMRZCDEewHgA6N69ewIvR/FSt4loTSfG0jYwUlHOtBF9UbqoGj6zFoJDrJt0XthY1KAwuE8nLN7gafHFIFTf777Aihcm+K9/c/l0rP35gBaP4dpmOB4iQGSduAOwlPJTIcQsAG8BOACgGkDY0S5SyjkA5gDKNqR4X4/io7WVJJTRwKO3Nuup84ZlSzlC6BZSxSp02jk0KAyauVr//UmJf7wyDYO/2gAA2NH+GJw7fg4aHOF/9bm2SUTJlFARlpTyOQDPAYAQ4kEAO8wYFJnHSDWwkcBTXulBpPMR+t77JuoPN6JLnguzxxRioklFWU6HwNThfSM+Ru8LxAm7vsIb/7jVf3396Hux+henIT/XiYO+Jh6QTkS2SrQK+hgp5S4hRHcAowD80pxhkVmiZbdGA0+0U5DUfbbqunBerjNq16looh1jqAqr1pUSc5Y+gAs/Xw8A2NUuH4Nufh4+h5JJqwGda5tEZKdEjyNc3LwG7ANwi5RyrwljIhPpbSUBjAc4ILb1Ua+vEQd9jXA6hOZ5t9HE2m1pcJ9OeHG90rP5+NqtWPX8H/z3jb9sClYd3/J7YcU3e3B/Cbs5EZG9Ep2CPsusgZA1SocWaK4B5+c6Y8r6IgVyLRKAr1EiRwBNEhGnrwGlkUajlDF9KQCUqfGXPtgOAHjy1VkYtvldAMAeV3sM/P0LONwqfNvSvOZgvWZzLTNgIrJNohkwpTg1qAT3RQaAvfW+iFuQQhvsh1YaG9UklYx29AA31myuhafOG3ew1Rrj5CUb0bN2G95+7mb/7b8fOQmv99FvIymhBGH1C0Gs27GIiMzAwxiyhN5pRurhBcH0mvCPHuBuEbhiEfw68Z6eE/p7Px1qwL2LZmJUzRoAwIHWLgz44zwcatU6jhFq/1kQESUi0mEMzICzRCx9ffX2+y6v/jau4Bv8OvE2Awn9vdZffYF1z9zov//W4aV47cRz4hxdyzEaxWP4iCgRDMBZIpa+vnqBKJFOV+rr6AX3aa/VRJwKDx77Q68/his2/hsAcDinFfpNWIBDzjaGx6K3Hh3LPuBEu4oREfE0pCwRS19fsxtSBL9OpOAe2pdaDXJq8O2+91tsnTXMH3xvv3Qiji8tjy34CmDcwO4J9zjmUX1ElCgG4CygZpFeX6P/nN38XCfatMrBxAVVYYcy6AXreAVvKYoU3EODV3CQu3/l/2HtnBv8950w8RUsOWlITONw5gjMvqIQ95f0w4xR/eDOc0FAWfuN9ZB5HtVHRIniFHSGC50qbZQSTofAgYMN/r7NodOnek344zly0J3nahHYIj1PaPDaWedF133f472//9Z/210X3YqF/S+MeRwAUHZ5f/9YEu1xzKP6iChRDMAZTmuqVKs5RuihDFoBKtYA7MwRKB1aEFas1K61w985K1ho8Hpo7XO4/P2l/usTJy5Cfev4AlzoF4FEae2vZjtLIooFA3CGi2VK1MzpUwHA1yQxfVlNWLbtzBFhXbJaBK/t24Hu3XF5831/GnoL5hde7H/eWCuxrQiMPKqPiBLFAJzhYulgFW361B3ludwaDTu0+kH7miTyXE60a9MqPHjdeSfwyCP+xy5f+yn+s24nRPPjYunGFTymspVbMHFBlamBkkf1EVEiGIAznF4rylBGssTSoQWYuKBKNwP96VCDvydzNPu8PlRNDVrL9XiUEuVmjw7/I544cSi6rNvZImDqNRQJlp/rxNThfVFS5OZ2ISJKWayCznAlRW6MHuD2Vz9ryXM5DVUBlxS5cUbvDrr3x7JPuEW2/ac/AV27+i9Pu3MR/nriUEgEAqZapa1VoR3MnedC5b0XtpgiTnS7UHmlB4NmrkavSSvCKsaJiOLFDDjDlVd6sHiDB40RWo7WeX3+gBStG9X/vk78wCt/tv3dd0DnzoE7Hn0Ugw71x66QDDe4QEwd3/RlNWHT21pZfKLbhZhBE5FVmAFnOK0MUEtopqll8pKP/cVUsRLN//j33L76TMvgu3s3MHFi1ICpVlTX1fuQn+tEnssZcS+v3rq20e1CbLhBRFZhBpzCYuk1XF7paZEV5rmcmDaib8zn+AZvRQp9fq+vKb43AqVyeevMS4Fdu4Bjjw3cMWsWcNdd/stI+2tDs9G99T64nA7MHlOo++eS6HYhNtwgIqswA05RwW0YtdZCQx9b+kp1iynZOq8PpYuqkds6tg5WeoEl0YzPIQTwwAMtg29tbYvgC0RumRlPNlpS5E6o61WiGTQRkR5mwCkqUrBRq3vV7Din+XzdUL4miQaNhheRBAeW0Kw6khwBtG/r1CzEyq/fh8onxgVuuP9+YMoUzeeJtL92osEOWlrPGe96LRtuEJFVGIBTVKSpT632knpiWbENDixqVq3VNUtLkwTatWmFdm1atZhCvvGDVzD5nbmBB37/PXDMMRGfSy9g2tH+kQ03iMgqQkb48DZbcXGxrKioSNrrpSOt4/eCuZuDjdGGFA6d7DhUnssJIYC6eh+65LlQf7jBUOar9Tw/HW5A7oEfUf3Xsf7bN98wAX3mzI75+YKFfvEAlC8NsR6kEPqcDK5EZBUhxAYpZbHWfVwDTiGhx++FUjNUowVAzhyBsad3i3qSUZ7LiUMNTdhb7/OvN8cTfAFl7fnaD5a2CL5vrPoo4eALJL6eGyqWdXYiIrMxAKeQSFuGgoONkSlXIYAxp3XzH72Xn+vUfJwzR0AIGNqqFE37gwewddYwTHn7WeWGSZMAKXHxBUUJP7eqpMiN0qEF6JLnws46L8pWbok7YHKLERHZiQE4hehltgLAuknn+TM9rUph9YADlZTA4g0elFd6UFLkRuW9F+KxMYUtAnGey4myy/ujLs5sN9hvNizDx49f6b8e+Pu5wIwZCT9vKDOzVm4xIiI7sQgrhRgtMtIqDNJaszVyxKD6PFqv63Lm4FBDEyL13jjiUD02PXaF/3rOqZfhwfN+61+rNlu06vBY6P155wiBXpNWcE2YiCzFAJwiyis9qD/cEHa73paX0GDaa9IKzec1UqxVOrQApYuqw7pcHW6UEYPvr6vewIMr/89/fcbNz2Nn+2Ms3aZjZtaqd1CFWrTGtpNEZCUGYAvEWlmrVd0LBLpZGfnw18vmRPPzR3sOrUrpRp3om3vYi09mX+6//seA4Zh+/o3+60QKo6IxcytS6EyC1n7qeLNrIqJouAZssnjWKPWKr9q1aWX4g790aAG0zjuSiN7FavqymoiZbrDLP17VIvieedNzLYJvfq7T0mCltf4tAAzu0ynq72qdalRS5Ma6Sefh65mXoklnuxbXhInICgzAJounstaMadWSIrdu041oz2Nky5Hr8EFsnTUMZW/8FQAwr/Ai9Lx7OXYcdWyLx53Y+UhD442Xerxi8JcNCeDF9dvQ9943db/oGPlixLaTRJRMnII2WTzB1KxpVZczR/PAhKNc2luQjLps02rMXvGo//rs8c9gW35nzceu/yrx4wojKa/04KUPtmt+2fjpcCNKX6kGAEPtOkOnl9l2koiSiQHYZPEEUzM++COdViS05qYNaOM7hI2PjUHrJqU4bGG/83HXJRMi/o6RrluxCA6iRzV32Yr0Gr5GiQkLqjAhpG+03u8EfzFi20kiSiYGYJPFE0yjffAbKeqKNMUdbZ9vfq4zbBp6+Cf/wRPLyvzXg294Gl93iB6IHHFGe633CKDFn6XWQQ+J0trixYBLRMnAAGyyeLMovQ/+0Appva0x0aa4I5k6vK//4IU2DYex4YlxOOKw8nxLTzwXE4ffGfH3g409vVvE+40EWvU9tmmVY0qHLj2cXiYiOzEAW8DMLMpo44lI25CMBBlnjsD5Ne/ib6/O9N825Ld/w5dHRw6oKodQ+k7fX9JP9zF6Xya0Aq3X12hJ8HUIgSYpOb1MRLZjAE5xRou6tKa+BYBxA7tH3YN878IN+O9j45B38AAAYHnBmfhDyaSoY8vPdWLqcGP7lAH9LxNWZrnBEj05ScUTlIjIDAzAKS6R9pRGAsP6R5/Dxy9O9V9feP2T+KxTz4i/444z6MS6n7Zdawd+OmxOcHYI0WI7WKInKEVbEiAiioYBOMXFUtQV09S3zwf06IGZ334LAFh13ECMv2xK1JJpd54L6yadZ/wNBNH7MqHHrOCbI8xrL2lmL2oiym5sxJHi4jkDV6vjUwvLlwOtWwPNwffi6/6K8aPuiRp8ja4n69HqYpUMoV2+EjlykCcoEZFZmAGngVgyW60p0tJF1Zi+rAb7DxzEe8/cgJ/V7QIAvH/8aRhb8mdDG4WNrCdrjSV0SnzGqH5he3TtEG/ANLMXNRFlN2bAGUZritTXJNFv03p8UTbSH3zXvPg6xl52r6Hg685zYfaYwogVzqH0Wj+mingDplYWz+1MRBSPhDJgIcREAL+D0o53I4DrpJQHzRgYRaeVYYZmdo6mRrz17M34+d6dAID3evTHVWPuBzZqd80KFe+ab6S1UodGW8hkSiRgslsWEZkl7gAshHADuBXAiVJKrxBiIYArAcw1aWwUgV417lEup79j1BlbqzB/wT3+3xl59SOo7hJb4Il3qjbSWum4gd3x4vptmvdrdeUyg4DyLTHeCu5g7JZFRGZIdA24FQCXEMIHIBfAzsSHREboZZhtnTlo5wCWzPk9Cn5Qgtz/up6IK349K66m0PFO1eqtleYIgXnrt8HlzMGhhiY0SWWL0MCf52Prbm9CxUx5Lieqpl6IQTNXh722GnzjreAmIjJb3GvAUkoPgIcBbAPwLYB9UspVZg2MItMLVAWbP0LNg8P8wXfUuDJcMe6huIJvIlO1ehXPjVJCAvD6mtCmlQNXDeyO9q5WWPflHv96cbz2NWf+rFQmonSQyBR0PoCRAHoBqAOwSAhxlZTyxZDHjQcwHgC6d++ewFApWF7IVK2QTVj2wkSc9P2XAIA9JxXizJEPoL4hvpAWa5crLcEtJoUAQpd9vb5GzFu/LaGgG0zN1lmpTETpIJEq6PMBfC2lrJVS+gAsAXBG6IOklHOklMVSyuJOnTol8HIULDiYFe+owdcPjfAHX6xdi+HjHok7+ALAQZ2jDY1Q16eDTy/Sq7kyK/gG71FmpTIRpYNE1oC3ARgohMgF4AUwBECFKaOiqP2G93l9ELIJi18sxSk7laYSnxzTC8OufRxfnXUWPCtWJPT6iXR30lqftlrwHuV0qVRmT2mi7BZ3AJZSfiCEeAXARwAaAFQCmGPWwLKZkX7DF/z4Neb87Y/+3xl75YN4v8fJcDdPs5qx1cfsCmir5LmcYXuUU71SmT2liSihKmgp5VQAU6M+kAxRMyKt9Ut/RlrYBTjnHMx5910AwOcdu2Ho9U+iKceBHADf7vOi56TEsl+V2RXQVpk2om/SXsss7ClNRGxFmQLKKz2Yvqwm6v7Xoz+tBnKG+K/XPTUfd+07FrLOC5czB15fk3mLqoi/77PWARJWyXM50zJgsVKbiBiAbRY6FalJSvxz4b04e2ulct27N7B5Mwa1aoV1zQ/pZVLWq8rPjT+wBa/BWpkJu5yOtMx+AVZqExF7QVsi6mlEQaa9VhMx+Pb9/ktsfWh4IPi+/jrwxRdAq5bfncxs7CgAXHpy54Seo6TIjXWTzsNVA83betautaPFGrc6ZRvpzzdVsVKbiJgBmyzSaUR19T5/tSugBN/grTotSInnX5mO875SCsvrf9YFudu2Ak5n2OvFe7SeHglg8QYPint0SHh6d8XH35ozKAD1hxvxwGUFGVG8lC6V2kRkHSGT2BS/uLhYVlRk9k4lrTaIoZw5AhCAr1H7z77Prq/x5j8CFc7rZz+PgROuC3tceaUHpYuq4Qs98NYkZrRuNKsgDIA/+9X682WbSSJKRUKIDVLKYq37OAVtMiNFNL4mqR18pcTflz7gD7617fJw3J1LcXXtsZrTrNNeq7Es+AKpVRCkNtpg8RIRZQoGYJPFW0RzXO032PrQcFz02fsAgBtL/oRT//AifA4nfI0S05fVhP2O7vS1SZJdEBSpW7WEMm2rNyYWLxFRumEANpneIQSR/PW1h/DW87cAAPa1aYfj71iKlQUtu3pacURfJGYVBDliOANCQj8Iq9PPLF4iokzBIiyThRbXtFX352rovXs73n72Zv/1LSPuxooTzjL8WjkCMGMGuk2rHMwafbLpBUH3lG+EzjK3LjUIB/9acIBl8RIRZQoGYAsEt0HUK8p6ZMWjGL1pNQDA26oNCm+dj0PONrrPmedyht0WLfheFeHg+2CHGposad340gfb4/o99exevQCb6m0miYiMYAC2WGjw7bnHg3eeudF/fduwO/Bq38FRn2dYf2VfbnAD/0j9ntX+yEYCsFXi7UXNimYiygYMwBYKrVye8cZfMfbjVQCABpGDkyYuxEFnW0PPtXiDx/9vdQ+sXoDLQaA/sttAX2at7DpWWif7xHMghNMhuJ5LRFmBAThBWn2c81xOTBvR198go1vdd3j36d/577/zkgl4pd/5Mb2O19eIlz7YrhnQgg+7V19bnaItHVqACQuqdJ/XmSMSbueod7LPzzvl4vNdP8X0XO1at+L0MhFlBQbgBJRXelD6SnXYnt46r8/fIOMvq/6G31QGmlGcMPEVeFsby3pD6WWTUiqFS3rrpZE6bvmapP+LQryBT+9kn1iDL6Ccc0xElA24DSkBZSu36HazOmbv99g6a5g/+N590R/R8+7lusE315kDAWXKWG9K2CH09/RIBDLP0KnvaSP6Rtwapfd7RpnZBIP7eYkoWzAAJ0Av8Pz57Wew7u/X+6/7TliIBf2Haj7WIQSuGtgdn9x3MWaPKQSgZNChodbldGDs6d2i7jFWDygIVlLkxoxR/eDOc0FAO5Br/Z5R8QRNp0MoLTmDcD8vEWUTTkEnIPRIuc4/1uL9vwV6Nv/5wt/jX0WXRHyOL2co94euowbvh3UHTS0X9+jgL3bSK2/S+mIQvHVH7+jCeDPZWM//dQcdSMH9vESUrRiAE1A6tAC3L6xCkwQmrXkeN/1vif++fhMWQLRvDxzSD0rBmajWOqoafIO35BjZYxwtIzX7LFp1PHcsrI5Y9SwAzB5TGLZGTUSUjRiAYxRa9XzM/t3431PX+O+fNmQ85haPUC4iBF8AGPjzfP/P8RwyoJV5GpnGjff3VOqWI0/QXmR3ngtjT+/WYptUMAFg3MDuDLhERM0YgGMQWvV859p/4g/vL/Tff/JtL+PHtkcYfr6tuwPBNZ6sNN62jIm0cwydKlczXk+dF4s3eDB6gBtrNteGBWdOLxMRtcQAHAO16rnTgb348P+u9t9+/+Dr8expo2J+vuDsNt6sNN62jPH+ntZUucrra8SazbWGulhpNe5ggCaibMIAbEDwlOtt783HxHXz/fcV3jofda72cT1vcHabjEMGzAh60Qq1jBRy6TXuALgmTETZgwE4CjVYuPbtwdYnxvlvn3XONfjbwMvjfl5nTnjLRSsPGTAr6OlNlQffH41e446ylVsYgIkoa3AfcBRlK7fgurUv4aOg4HvKH+clFHyFAMac1i2pwSZS0ItFpPOOjRZyxVNwRkSUaZgBR7J7N9ZNHuK/fOTMcXhi0FjNh141sDvuL+kHAOips89WJaVyqEJxjw5JC8JmBb3gqfJ4C63M3gZFRJSOGID1PPwwUFrqvyz+w7/wQ7t83Ycvr/7WH4CNnAKUjCnX4DXfHJ0xxRP0Ep0qT3QbFBFRJmAADrV3L9Chg/9yy29vRUnni6N2earz+lBe6UFJkdvwEXyJTrneU77Rf0KSQwiMPb2b/0uA3nahYHYFvWQUnBERpToh4zw0PR7FxcWyoqIiaa8Xs8cfByZMCFzv3Al07mwokwQCXasKp6/SPX0omEMIPHJF/5gDjxJcP4bX1xR2X5tWOTjUEH578Gs2ScmgR0SUBEKIHmmjRgAACbVJREFUDVLKYq37mAEDwL59QF5e4Pquu4BZs/yXwVOu5ZUe3fN1d9Z5UV7pwU+HGwy9bKOUKF1U7X8NIwKZrXaQjRR8AaBJSnw981JDr0VERNZhBvzUU8AttwSut28HunaNuGe26C+r/K0og7mb11MjbdPR4nLmoEO7NoamY/X6PxsV2luaiIisEykDztptSMvf26LsB2oOvl+Mu0EpT24OvqWvVMPTfOKQp86L0leq/eflTh0efr6uup4az7qu19fU4rUmLKjCiX9+A0V/WYVek1Zg0MzV/tdOZN2YhU5ERKkjKwNw5T2zMOysPv7rM25+HsN7jfIHuenLavz9nlW+Ronpy2oAhJ+v685zYcaofigpcpu2labe14S99T5/UJ68ZCPKKz1xP3/wGImIyH7ZNQV94ABw5JH+yxdOuRRTL7i5xUOibSHaGmX9NLT6GFAyT6Nn5Uai7rXVO3s3RwBNIUN35giUXR57oRcRESWOU9AAMHdui+B75o3PhgVfQHu7TjA1S9ajlx0Hn/0br511Xs3nf2xMIbbOvBSPXlGI/Fyn//F5LieDLxFRisr8DLi+HjjiCGV9FwBuuAGDfn5l3IVM8RQxlVd6MGXpRvx0OLEsmAVURETpJaO3IUU84WfePOCqqwIP/vxz4Be/QKnGNLFRnjovBs1cbbiBhNaUdDxYQEVElFnSOgDrnfCTc9CLEYNPAg4dAgC8fsqFuOWCW9HllW0Y3MeLNZtr4fU1QohAYmyUQGCbUaQThcorPZi+rEZzu5IReS4nhADq6n1smkFElIHSOgBrnfAz5OM1GHH/xf7ri29+Bp+27wxACZgvrt/mvy+e4Bv6K1o9ndVtTKGV1Ea4nA5WKxMRZYG0DsDBe2JbN/iw4Ylf48jDzbeNHYtBJ/8uoaYVwSJly6F7c8tWbokr+AJg8CUiyhJxV0ELIQqEEFVB//wohJgQ/TfNo+6JvWjLOnz2yGX+4DtuwnPA/PmmnS+bg8jZcuje3Hhf1yEEgy8RUZaIOwOWUm4BUAgAQggHAA+ApSaNy5C7zuuFs88tRL73RwDA68efgTuu+DNmjFJOBNI7dzZWkborhxZHlVd6Ih7YEMnY07vFMToiIkpHZk1BDwHwpZTyG5OeL6LySg/mvbQGi8qu9t920fVPYv9xJ2BGULFS6dCCuNdijQg9hF4tCosn+F41sLv/KEEiIsp8ZgXgKwG8ZNJzRaQGucGfKW0h3/rFabh1zDTMGH2y9vStRductfbkahWFqfJcTvgamzT3AufnOhl8iYiyTMIBWAjRGsAIAJN17h8PYDwAdO/ePdGX8we51/uciZ59lis3NjSFVSKrj/WF9maMkSNHIAdo8Tx6e3L11n4FgKqpF2pWRzsdAlOH901ojERElH7MaEV5MYCPpJTfa90ppZwjpSyWUhZ36tQp4RfTC3Jat5tRhPXI5f1Rdnl/zYMXQukdlKDeXlLkRtmvWj5X2a/YKpKIKBuZMQU9Fkmafgb0C6u0gl+iRVjuPJc/OBoJkloHJYRmyyVFbgZcIiJKLAMWQuQCuADAEnOGE13p0ALds3hDDe7TCfEegeB0iJhbP0Y6ppCIiChYQhmwlLIeQEeTxmKIGsx0+z83K6/0YPEGT1w1WPm5Tkwd3jeuwMkMl4iIjEjLTlhGglykimQtDiHwyBVcjyUiouTI2POAYynAcjkdDL5ERJRUGRuA9SqSVQ4huE5LRES2ScspaCNKhxZgwoIq3fuZ8RIRkZ0yNgMuKXIjP9epeV+ey8ngS0REtsrYAAwAU4f31dyyNG0EO08REZG9MnYKGjC+ZYmIiCjZMjoAA9yXS0REqSmjp6CJiIhSFQMwERGRDRiAiYiIbMAATEREZAMGYCIiIhswABMREdmAAZiIiMgGDMBEREQ2YAAmIiKyAQMwERGRDYSUMnkvJkQtgG+S9oKxORrAD3YPIsn4nrMD33P2yMb3nervuYeUspPWHUkNwKlMCFEhpSy2exzJxPecHfies0c2vu90fs+cgiYiIrIBAzAREZENGIAD5tg9ABvwPWcHvufskY3vO23fM9eAiYiIbMAMmIiIyAZZH4CFEBOFEDVCiE1CiJeEEG3tHpPVhBC3Nb/fGiHEBLvHYxUhxPNCiF1CiE1Bt3UQQrwlhPi8+d/5do7RbDrv+fLm/9ZNQoi0rBaNROc9lwkhNgshPhZCLBVC5Nk5RrPpvOf7mt9vlRBilRCii51jtILW+w66704hhBRCHG3H2OKR1QFYCOEGcCuAYinlSQAcAK60d1TWEkKcBOAGAKcB6A9gmBDiOHtHZZm5AC4KuW0SgLellMcBeLv5OpPMRfh73gRgFIC1SR9NcsxF+Ht+C8BJUsqTAXwGYHKyB2WxuQh/z2VSypOllIUAlgO4N+mjst5chL9vCCG6AbgAwLZkDygRWR2Am7UC4BJCtAKQC2CnzeOx2gkA1ksp66WUDQD+A+Aym8dkCSnlWgB7Qm4eCeCF5p9fAFCS1EFZTOs9Syk/lVJusWlIltN5z6ua/34DwHoAXZM+MAvpvOcfgy7bAci4Ah+d/6cBYDaAu5Bm7zmrA7CU0gPgYSjfmr4FsE9KucreUVluE4CzhRAdhRC5AC4B0M3mMSXTsVLKbwGg+d/H2Dwest71AN6wexDJIIR4QAixHcA4ZGYGHEYIMQKAR0pZbfdYYpXVAbh5/W8kgF4AugBoJ4S4yt5RWUtK+SmAWVCm6N4EUA2gIeIvEaUpIcQUKH+/59k9lmSQUk6RUnaD8n7/YPd4rNacRExBmn7ZyOoADOB8AF9LKWullD4ASwCcYfOYLCelfE5KeYqU8mwo0zmf2z2mJPpeCNEZAJr/vcvm8ZBFhBDXABgGYJzMvv2W8wGMtnsQSdAbSgJVLYTYCmWp4SMhxM9sHZVB2R6AtwEYKITIFUIIAEMAfGrzmCwnhDim+d/doRTnvGTviJLqNQDXNP98DYBXbRwLWUQIcRGAuwGMkFLW2z2eZAgpphwBYLNdY0kWKeVGKeUxUsqeUsqeAHYAOEVK+Z3NQzMk6xtxCCGmAxgDZZqqEsDvpJSH7B2VtYQQ7wLoCMAH4HYp5ds2D8kSQoiXAJwL5bSU7wFMBVAOYCGA7lC+gF0updQq6khLOu95D4AnAHQCUAegSko51K4xmk3nPU8G0AbA7uaHrZdS3mTLAC2g854vAVAAoAnKqXM3Nde5ZAyt9y2lfC7o/q1QdrWk8ulIflkfgImIiOyQ7VPQREREtmAAJiIisgEDMBERkQ0YgImIiGzAAExERGQDBmAiIiIbMAATERHZgAGYiIjIBv8Plmtjp4IwxQwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     }
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "source": [
    "df = pd.DataFrame({'Stationary Series':y-0.8602*x, 'Mean':np.mean(y-0.8602*x)})"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "source": [
    "def zscore(series):\r\n",
    "    return (series - series.mean()) / np.std(series)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "source": [
    "zscore_calcu = zscore(y-0.8602*x)\r\n",
    "print(zscore_calcu)"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "date\n",
      "2015-01-05    1.746636\n",
      "2015-01-06    1.577667\n",
      "2015-01-07    1.828576\n",
      "2015-01-08    1.533107\n",
      "2015-01-09    1.721811\n",
      "                ...   \n",
      "2017-07-12   -1.025160\n",
      "2017-07-13   -1.000172\n",
      "2017-07-14   -0.838614\n",
      "2017-07-17   -0.375674\n",
      "2017-07-18   -0.482295\n",
      "Length: 619, dtype: float64\n"
     ]
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "33347927775e93bf1428669ae9a55f1131bd6123ff01e11559ab145ac8d0794b"
  },
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.10 64-bit ('user': conda)"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}